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Abstract. On large scales a nonlinear transformation of matter density field can be viewed as 
a biased tracer of the density field itself. A nonlinear transformation also modifies the redshift 
space distortions in the same limit, giving rise to a velocity bias. In models with primordial 
nongaussianity a nonlinear transformation generates a scale dependent bias on large scales. 
We derive analytic expressions for the large scale bias, the velocity bias and the redshift 
space distortion (RSD) parameter f3, as well as the scale dependent bias from primordial 
nongaussianity for a general nonlinear transformation. These biases can be expressed entirely 
in terms of the one point distribution function (PDF) of the final field and the parameters 
of the transformation. The analysis shows that one can view the large scale bias different 
from unity and primordial nongaussianity bias as a consequence of converting higher order 
correlations in density into 2-point correlations of its nonlinear transform. Our analysis 
allows one to devise nonlinear transformations with nearly arbitrary bias properties, which 
can be used to increase the signal in the large scale clustering limit. We apply the results 
to the ionizing equilibrium model of Lyman-a forest, in which Lyman-a flux F is related to 
the density perturbation 5 via a nonlinear transformation. Velocity bias can be expressed 
as an average over the Lyman-a flux PDF. At z = 2.4 we predict the velocity bias of -0.1, 
compared to the observed value of —0.13 ± 0.03. Bias and primordial nongaussianity bias 
depend on the parameters of the transformation. Measurements of bias can thus be used to 
constrain these parameters, and for reasonable values of the ionizing background intensity we 
can match the predictions to observations. Matching to the observed values we predict the 
ratio of primordial nongaussianity bias to bias to have the opposite sign and lower magnitude 
than the corresponding values for the highly biased galaxies, but this depends on the model 
parameters and can also vanish or change the sign. 
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1 Introduction 

A simple model for Ly-a forest relates the neutral hydrogen responsible for absorption to the 
underlying matter density field via a chemical equilibrium equation, where recombinations 
and photo-ionizations balance each other [1]. Optical depth is proportional to neutral hydro- 
gen density, which in the ionizing equilibrium can be related to gas density as r = A(l + 5) a 
. Here 5 is the gas overdensity parameter, a = 2 — 0.7(7 — 1); where 7 — 1 = din p/dlnT 
is the slope of temperature-density relation and the recombination coefficient is assumed to 
scale as T~ ' 7 . Typical value is a = 1.6, with a = 2 being the isothermal case. The observed 
flux F is related to the optical depth r as F = exp(— r). The relation between the observed 
flux and the density perturbation 5 is thus highly noninear. 

Even though the relation between the Ly-a forest and the underlying matter density 
is nonlinear, on large scales Ly-a forest fluctuations trace the dark matter fluctuations up 
to a constant factor of proportionality called density bias. The basic premise of Lyman-a 
forest clustering analyses is that the bias is a known function of the underlying parameters. 
This is the basis of the statements that Ly-a forest measures directly the amplitude of matter 
fluctuations at 2 < z < 4 [2] . However, so far all of the predictions came from simulations and 
we do not have a good analytic understanding of how the bias in Ly-a forest is determined 
by the parameters of the model. The purpose of this paper is to derive the bias analytically 
and to explore its sensitivity to the physical parameters. We will assume the gas density is 
related to the dark matter density smoothed on the Jeans scale (or, more precisely, filtering 
length, [3]). This relation is not exact and there is scatter around it, at the level of 10-50% [4]. 
Initially we assume the relation between optical depth and matter density is deterministic, 
later we generalize this to include a simple form of scatter. 

1 In some previous papers P is used in place of a, but here we will reserve /3 for the redshift space distortion 
parameter. 
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Our second motivation is to derive an analytic prediction for the redshift space dis- 
tortions (RSD). In redshift space the observed position is a sum of the radial distance (in 
velocity units) and the radial velocity. The velocity gradients give rise to additional per- 
turbations called RSD. One can view this as a mapping from the real space to the redshift 
space with the total number of tracers being conserved. In the large scale limit the RSD take 
a particularly simple form first derived by Kaiser [5]. In this limit the large scale velocity 
of any tracer follows the dark matter and there is no velocity bias. However, if the field 
is transformed after the RSD mapping then the transformed field acquires a velocity bias 
different from unity. This is the case for Ly-a forest since RSD act on optical depth r, while 
the observable is F = exp(— r). In this paper we derive its velocity bias. Combining with 
the bias predictions one can also predict the RSD parameter /3, which is defined as the ratio 
between the velocity bias and the density bias. Recent observations of bias and f3 provide an 
opportunity to compare our predictions to the observations [6] . 

The third motivation for this paper is to explore the sensitivity of Ly-a forest, and 
nonlinear transforms in general, to the primordial nongaussianity. Primordial nongaussianity 
models have an additional contribution added to the primordial density field, which does not 
show up in its 2-point correlation function, but only in the higher order correlations. However, 
such a component can show up in the 2-point correlations of the nonlinear transform of the 
density field. We derive its amplitude and scale dependence and show that it agrees with the 
corresponding scale dependent bias of biased halos [7-9] up to a prefactor which is determined 
by the transformation parameters. 

Although the Ly-a forest will be our primary application in this paper, the formalism 
we develop here is more general and can be applied to any tracer of the dark matter. The 
goal of this work is to investigate the large scale density bias (hereafter bias), velocity bias 
(as defined in redshift space distortions) and scale dependent primordial nongaussianity bias 
of a general nonlinear transformation of the matter density field 5. For example, galaxies are 
formed inside dark matter halos and these are often modeled as a nonlinear transform of the 
local density field, 5h = bid + &2^ 2 + ••• [10], where density is smoothed on a scale typically 
related to the Lagrangian scale of galaxy sized halos (which is of order 1 Megaparsec). The 
linear bias b\ is a monotonically increasing function of mass of the halos in which the galaxies 
live. This is usually explained within the context of a universal halo mass function and the 
peak background split [11]. One can however also view the bias different from unity as a 
consequence of a nonlinear transformation of the density field: as we will show in this paper a 
nonlinear transformation can result in an arbitrary value of bias. This analogy is not perfect: 
galaxies and halos are discrete objects and in addition to the nonlinear transformation one 
must also impose the exclusion constraint, where no halos can be within the virial radius 
of each other. However, one can also consider further nonlinear transforms of the galaxy 
field which modify the bias properties. This may be useful for observations where individual 
galaxies are not measured, only their overall intensity imprints, such as in the 21-cm intensity 
mapping [12]. 

2 Formalism 

We begin with the derivation of the large scale bias. Let us call the nonlinear transformation 
of the density field t(5) and decompose the density perturbation into a long wavelength 
component 5i and a short wavelength component 5 S , 5 = Si + S s , both with zero average. We 
will also assume \6i\ <C \5 S \ in an rms sense. This assumption is well justified in our universe 
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where power per mode k 3 P(k) is rapidly increasing with wavevector k (P(k) is the power 
spectrum) . 

We want to know the response of the nonlinear transform r(<5) to a long wavelength 
mode 61. We define the (density) bias as 

Here () denotes average over the field. We can expand t(5) in a Taylor series 



ni 

n=0 

where r^ n '(0) is the n-th derivative of function r evaluated at 5 = 0. The first term in 
expansion above is a constant and does not depend on 5i. The second term, n = 1, is linear 
in (5/, so its dependence on 5i is simply 7"W(0)(^ and the first order bias of r is b\ = r^\0). 

To understand the bias for higher order terms (n > 1), we need to develop understanding 
of gravitational coupling between the long wavelength and the short wavelength modes. This 
can be achieved using perturbation theory, but here we will pursue a simpler approach in 
terms of a constant overdensity 5i. A similar derivation has recently been given in [13]. 

2.1 Coupling between long and short wavelength modes 

We would like to analyze the response of small scale perturbations to a long wavelength 
density perturbation 5i, which we will model as constant in space. This leads to a locally 
slightly overdense or underdense universe, but one where the global time and coordinates 
are still given by the global unperturbed value. Thus one must consider a slightly changed 
Hubble expansion rate, as well as the fact that the small scale density perturbation feels the 
additional gravity from the overdensity (or underdensity) of the long wavelength mode. In 
addition, the long wavelength overdensity itself grows in time according to the linear growth 
rate. Most of the time during matter domination the universe is EdS with Q m = 1, so we 
focus on that solution. We write the expressions in terms of an overdensity perturbation 
5i > 0, although the final result is the same for an underdensity <5/ < 0. 

We begin with the standard approach in deriving the spherical collapse model solution. 
We work in the locally comoving coordinates, but using global time t. In a homogeneous 
universe the cycloid solution to the expansion is 

— = - (1 - cos rj) , = - (77 - sin i]) . (2.3) 

Taylor expanding in small r/ and eliminating r/ from the expressions order by order consistently 
to get the first two non- vanishing terms gives 

£) 2/3 (l-i^ i0 ). (2.4) 

Here time to is the age of the EdS universe today, and the factors have been arranged such 
that the long wavelength mode is 5i = 5(a~ 3 ) = t 2 ^ 3 5io, i- e - the long wavelength mode is 
growing according to its linear growth rate, which in an EdS universe is just the expansion 
rate a = (t/to) 2 ^ 3 . We have applied mass conservation, so the long wavelength density 
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perturbation is given simply by the change in volume 6 (a 3 ). The corresponding Hubble 
rate to the same order is 

HH'N- < 2 - s > 

A short wavelength perturbation 5 S in this universe obeys the equation 

5s + 2% = 4irGp5 s (l + i 2 / 3 ^) . (2.6) 

This is the usual second order equation derived from continuity and Euler's equation for dark 
matter. The last term includes the fact that gravity responds to the total matter and so in 
a slightly overdense region (caused by 5[) the gravitational force will be slightly stronger. 
Note that this term vanishes initially, at t = 0, so it does not enter in setting up the initial 
conditions and there is no coupling between the modes initially, as expected. Here p is the 
global density and we have H 2 = 8nGp/3 = (2/3t) 2 . We know that for 5o = the solution 
is 5 S = 5 s ot 2 / s , so we can write the ansatz solution as 5 S = 5 S Qt 2 ^(l + (3 2 t 2 ^ 3 5io). Inserting 
this ansatz into equation 2.6 gives f3 2 = 13/21. 

This solution was with respect to the local expansion, i.e. wrt to the local comoving 
coordinates. With respect to the global coordinates we have 

S 9 S = (1 + 605s = (1 + 5^(1 + p 2 5i)a5 s0 ~ f 1 + ^SiJ a5 s0 = (1 + u 2 6 t ) a5 s0 . (2.7) 

We thus derived the result that the short scale density perturbation is enhanced in the 
presence of a long wavelength perturbation 5[ by a multiplicative factor proportional to 5[. 
The coefficient of proportionality v 2 = 34/21 is the well known angular average of the second 
order perturbation theory kernel ^(kj^ka) [14]. The small scale perturbations are also 
rescaled in size by 5i/3, but this effect is not relevant for the purpose of this paper. 

2.2 Large scale bias 

With this result in hand we can write to the lowest order in 6i, 

5 n = (5 S + S t ) n = [6.(1 + v 2 6i) + 5i} n ~ 6^(1 + nv 2 5{) + n6^ x 5 l . (2.8) 

Thus 

3r ^T^(0)d6 n ^ T^(0)5 n ^ T^(0)5 n - 1 I Ar\ I dr 
b T=dsr^^^^ = U2 ^ n ^^ + ^ n n\ =U2 \ S d6/ + \d6 

n=l n= 2 n=l 

(2.9) 

We used 5 ~ 5 S at the lowest order. 

The simplest non-trivial example is that of a quadratic dependence, such as the optical 
depth of Lyman-a forest in the case of isothermal density-temperature relation, r = A(l+5) 2 . 
We find b T = 2A(1 + v 2 (T 2 )^ where aj is the rms density field smoothed on a Jeans scale. 
Note that we can get the variance from the transformed field itself, i.e. a 2 = A~ 1 (t) — 1. 
A more meaningful way to express this is to look at the bias of optical depth overdensity 
5t = t/(t) - 1, 



2(1 + 



= i, „2 • ( 2 - 10 ) 
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In the limit of small a 2 this agrees with the linearized analysis, b$ T = 2, where only the 
first term in Taylor expansion is kept. On the other hand, in the limit of large a 2 we get 
bs T = 2u2 = 3.24. In the case of Lyman-a forest applications, for 2 < z < 4, a 2 is expected 
to be of order unity, which means that the bias contribution from the quadratic term is 
important relative to the linear term. The derivation assumed that the short wavelength 
modes grow according to linear theory, so this approximation presumably breaks down if the 
variance of the density field aj S> 1. Limited tests in simulations however suggest that this 
approximation works well even if aj ^> 1, e.g. figure 1 of [15]. 

2.3 Primordial nongaussianity 

The case of primordial nongaussianity of local type is even simpler. The local model for 
initial potential is 

$ nG (x) = <p(x) + / NL (^ 2 (x) - (if 2 )) , (2.11) 

where (p is primordial Gaussian potential. Performing the same long-short wavelength mode 
split as above one finds, 

x X (1 L o a „ , 

hi 

where the relation between initial potential and final linear density is given by 



S,*G = Ss(l + 2f mj a£ r 6 l ), (2.12) 



2k 2 c 2 D(z)T(k) 

"/nl(M) = ^rr2^ > ( 2 - 13 ) 

O-flg " 

where Hq is the Hubble parameter, T{k) the transfer function an D{z) the growth rate. Note 
that «/ NL scales as k 2 on large scales where the transfer function is unity, hence the effect 
becomes large on very large scales. 

We see that in equation 2.12 replaces v-i in equation 2.7. The large scale bias 

due to the nongaussian term is thus 

KNG = 2h L aJ^(6^y (2.14) 



NL 



Other primordial nongaussianity model can be treated the same way, except that a/, 
changes, e.g. it scales as k" 1 for orthonormal models and is constant for equilateral models 
[16]. 

2.4 Redshift space distortions 

A third application are the redshift space distortions(RSD). The situation we wish to con- 
sider is one where observations are in redshift space, meaning the position of the object (or 
intensity for continous case) has to include the peculiar velocity of the object. There are two 
cases that can be considered. One is when RSD transformation occurs after the nonlinear 
transformation. In this case RSD simply remaps the variable r from the real space to the 
redshift space, such that the total is conserved, i.e. r(r)d 3 r = r(s)(i 3 s, where r denotes real 
space coordinate and s the redshift space counterpart. The Jacobian of the transformation 
is | d 3 s / ti 3 r | = 1 — dv z /dz, where we denote with z the radial direction and v z is the velocity 
in the radial direction. If we consider the response to a Fourier mode 5i then the usual linear 
order result gives dv z /dz = ffj?5i, where \i is the angle between the Fourier mode direction 
and line of sight and / is the logarithmic growth rate [5]. At linear order we thus have 

r(s)=r(r)(l + / / x 2 ^). (2.15) 
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we thus find that the lowest order contribution from RSD is (r)/^ 2 <5;. Note that if the field 
is divided by the mean then there is no velocity bias in RSD, i.e. velocities are a faithful 
tracer of the matter field, which is the usual result for galaxies. 

A more complicated case is where the nonlinear transformation occurs after RSD map- 
ping. An example is the flux in Lyman-a forest. The optical depth r is proportional to the 
neutral hydrogen along the line of sight, which is remapped into the redshift space due to 
peculiar velocities, giving optical depth r(s). The observable however is the flux F relative to 
unabsorbed value (continuum), i.e. the fraction of the flux absorbed by the neutral hydrogen 
is -F[t(s)] = exp[— t(s)]. To preserve generality we will however assume i ? (r)(s) is a general 
function of r. We can expand in r, and at the lowest order in 5i we have 



oo 



F[r(s)] = F[r(r)(l + //A)] = ^ v ' v , ^ " = F[t(t)] + fn% (r 

n=0 ' ^ 

(2.16) 

Combining all previous results we thus arrive at the final expression for the bias, 
t , , , 8F[t(S(b))] ( n , i\/ r dF\ /dF\ , ,/ dF\ fnirr . 

This is the central equation of this paper. It has a remarkable property that if the nonlin- 
ear transform is monotonic we can invert the relations and express everything in terms of 
one point distribution function (PDF) of the final observable F and the parameters of the 
nonlinear transformation. Often the fluctuating field is normalized to unity, i.e. the field is 
defined to be 8F = F(s)/{F) — 1. In this case the bias above is divided by (F). 
To simplify the notation we can define the (density) bias of F as 



the primordial nongaussianity bias as 



and velocity bias of F as 



h,-(r^) (2.20) 



such that 

Mm, /nl) = b F + b fNh f NL aJ^ + fK/j 2 . (2.21) 

It is customary to introduce the RSD parameter /3 as 6f(/nl = 0) = bp (/J- = 0)(1 + /3/i 2 ), in 
which case ^ 

°F M d ds) + ids) 

2.5 Galaxies and nonlinear transformations 

We have shown above that a nonlinear transformation of the density field changes the bias. 
This happens even if there is no bias at the linear order, i.e. in 5h = b±S + bi& 2 ■■■ the 
unrenormalized value is b\ = 1. It is well known that dark matter halos in which galaxies 
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form have bias different from unity, ranging from 0.7 at the low mass end to an arbitrarily 
large value at the high mass end [11]. In the picture pursued in this paper galaxies and 
halos with b 7^ 1 can be viewed as a result of a nonlinear transformation. In this picture the 
bias different from unity happens because higher order correlations in 5 show up as 2-point 
correlations in the nonlinear transform of S, hence the large scale bias is renormalized [17, 18]. 
So one can ask how far can this analogy be pursued. For example, we have seen that for 
t = 5 n one gets b$ T = nv-z if n is even, so one can obtain an arbitrarily large bias from such 
transforms. One can lower this to an arbitrary low number, for example by adding a constant 
component, e.g. for r = A + 5 n one has 

bsr = -t^tt, (2.23) 

for even n. As we have seen in previous section in the application to Ly-a forest one can 
also easily obtain |6| < 1. 

Primordial nongaussianity also leads to a scale dependent bias: for r = 5 n we find 
6/ NL = 2n, so we have 

1.2. (2.24) 



If there is a linear component added to this with a positive sign, e.g. r = 5 n + AS and 
A > 0, then this number will be reduced, because a linear component adds to the bias b$ T 
but not to the primordial nongaussianity bias. Of course, if only the linear component is 
present then bf NL = 0. The ratio bf NL /b can also be increased if we consider A < 0, which 

reduces b$ T but not 6/ NL . As discussed above this can be compared to = 25 c (b — l)/b ~ 
3.37(6 — l)/b, where S c = 1.68 for biased halos [7]. We thus find that a nonlinear transform 
can accommodate the complete range of allowed values of bias and primordial nongaussianity 
bias. 

A few more remarks are in order. First, if these nonlinear transforms act on the den- 
sity field and the redshift space distortions act on the nonlinear transform then there is no 
velocity bias, consistent with the standard picture of halos and galaxies. Second, halos are 
discrete objects and in addition to the nonlinear transformation one must also impose the 
exclusion constraint, where no halos can be within the virial radius of each other. This con- 
straint cannot be expressed in terms of just the local density transformation (the situation 
is even more complicated for real galaxies, which can be either central galaxies or satellites 
in a halo). Third, even if the nonlinear transform can give an arbitrary value of the bias, 
the corresponding shot noise can differ significantly from that of the discrete halo tracers. 
So in detail there are considerable differences between the real galaxies or halos and the 
nonlinear transform of the density field, even if both display large scale bias and primordial 
nongaussianity bias. 

One can also consider further nonlinear transforms of the galaxy density field, which 
will result in a bias different from the original field. For example, this way one can construct 
a tracer with a different large scale bias than the original galaxy field, which may be useful 
if one wants to apply the multiple tracer sampling variance methods as in [19, 20]. One such 
application where this could be useful is 21-cm intensity mapping [12]. 

3 Application to Lyman-a forest 

The above derivations are entirely general and can in principle be applied to any nonlinearly 
transformed density field. Let us apply the above results to the model for Lyman-a forest, 
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ignoring the primordial nongaussianity for the moment. The optical depth is given by r = 
^4(1 + S) a and the flux is F = exp(— r). This gives 

b F = -Aau 2 <(1 + 5) a - l 5F) - Aa <(1 + S^F) . (3.1) 

This can be rewritten as 

b F = -Aa(u 2 - 1) <(1 + 5) a ~ 1 5F) - Aa ((1 + 5) a F) . (3.2) 

The velocity bias is 

b v = -A ((1 + 5) a F) , (3.3) 
The redshift space distortion parameter 13 f is 

f^ = a + {u 2 -l)^L. (3.4) 

Note that for z > 2 we can approximate / ~ 1 to a high accuracy and we will assume 
/ = 1 below. In the expression for /3p above, the term ((1 + 5) a F) is always positive, while 
(5(1 + 5) a ~ 1 5F'j can be of either sign because 8 is of either sign. As shown below it is 
typically negative, and as a result the RSD parameter is typically larger than a" 1 . 

Since the relation between matter density and flux is monotonic we can invert the above 
relation, i.e. we can write 5(F) = (— In F/A) a — 1. This way we obtain an expression that 
only depends on the parameters of the nonlinear transformation, A and a, and on the one 
point distribution function (PDF) of the final field F itself. In terms of flux and its PDF 
bias can be expressed as 

bp = a{FhxF) + a{u 2 - 1) ^FlnF[l - (- lnF/A)-^ 1 }"} , (3.5) 

and 

b v = (FlnF). (3.6) 

We can thus use the physical model parameters and the observed flux PDF to determine the 
bias, without any need for simulations. Note that velocity bias b v is completely independent 
of the parameters of the nonlinear transformation and is negative definite across the entire 
range of F, hence it is always negative. In contrast, bp explicitly depends on A and a, so 
these two parameters must be specified in addition to flux PDF. Moreover, the last term 
contains (F(l + 5) a ~ 1 ), which becomes comparably large in the voids where 5 ~ —1, and this 
terms contributes with opposite sign to the first two terms in equation 3.5. The predictions 
for bp are thus very sensitive to the void regions. 

3.1 Log-normal models predictions 

To proceed we must evaluate these terms. We could simply use the observed PDF of F 
and explore the predictions as a function of A and a and we will do so in the next section. 
However, we also want to develop a better analytic understanding of where the dominant 
contributions come from, so we will use the log-normal model, which we show to give a 
reasonable approximation to the PDF of the flux field. In this model the nonlinear density 
field is given by 1 + 5 = exp(<5c — cr 2 /2), where 5q is a gaussian random field, with a gaussian 

probability distribution p(5c|cr) = (27TO" 2 ) -1 / 2 exp[— p], Note that the variance a is related 
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to the Jeans smoothed variance of density field, a 2 = exp(cr 2 ) — 1. For a given value of a 
the constants A and a can be determined by normalizing to the observed mean flux and rms 
variance of flux. The mean flux is F = 0.47 ± 0.02 at z = 4, F = 0.68 ± 0.02 at z = 3 and 
F = 0.82 ± 0.01 at z = 2.4, while ap = 0.013 ± 0.003 at z = 4, 0.0117 ± 0.006 at z = 3 and 
0.0079 ±0.007 at z = 2.4 [21]. 

With this the model is fully specified and the bias and (3 can be calculated by doing 
two simple gaussian integrals. These can be written as 



e e z^aoc (3.7) 



and 



S) a , 



and 



bv = ((1 + 5) a F) = (2TTCT 2 )- 1 / 2 J e «(5c-«y2) e -Ae^G-° /*) ^ ^ g) 

These integrals can be written in a simpler form by transformation dF/dSc = —AaF(l + 



b F 



dFp(5G\cf) 



1 + 5 



b v = —a 1 / dFp{5c\(y) 
Jo 



(3.9) 



(3.10) 



and we note that 5g is an implicit function of F whose explicit form can be obtained by 
inverting the relation F = exp[— Ae a ^ G ~ a Z 2 )]. 

Finally, for a given a the parameters A and a can be determined from the mean flux 
and rms variance, which are given by 



and 



(F) 
(F 2 ) 



1 



dFp(5 G \a)- 
Aar 



dFp ^A^r' 



with a 2 F = (F 2 ) — (F) 2 , or 

a 2 F = C dFp(5 G \a) 
Jo 

We can write these equations above as 



(F - (F)f 
AarF 



X = J dFWx(F), 



(3.11) 
(3.12) 

(3.13) 

(3.14) 



where Wx is the integration kernel and X = bp, b v , F, a 2 ? and can be either taken from the 
lognormal model or directly expressed from the PDF, as in equations 3.5-3.6. These have 
a simple interpretation: the measured value of X is simply the average of the integration 
kernel over the interval < jP < 1. 

For z = 2.4, we find that a = 1.5, A = 0.3 and a = 1.6, give a reasonably good fit to 
the flux PDF, as shown in the bottom of figure 1. The corresponding mean flux and rms 
are (F) = 0.82 and a 2 ? = 0.07, in good agreement with observed values [21]. The predicted 
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Figure 1. Top panel shows window functions Wx for bias (Wb F ) and velocity bias (b v jr) of F, 
as well as for mean flux (Wf) and flux rms (W^a ) for a representative log-normal model (a = f.5, 
a = 1.6) at z = 2.4. The window has the same sign for velocity bias b v p over the entire interval, while 
the window for bias Wb F crosses zero and becomes large for F > 0.93, suggesting a large contribution 
from void regions. We also show window functions of H = F — F 2 /2, both bias (Wb H ) and velocity 
bias (b v n), for which the void contribution is suppressed. Bottom figure shows the PDF from this 
model, compared to observations in [22, 23], indicating that log-normal model provides a decent fit to 
the data. We also show the deconvolved PDF of dn/dF derived from the log-normal PDF of dn/dF, 
assuming a scatter of a e — 0.4, as described in the text. We see that this amount of scatter has a 
relatively modest effect on the PDF. 



values of bias are b v = —0.1 and bp = —0.11. At higher redshifts the agreement between 
observed flux PDF and log-normal model becomes worse, so log-normal model is less useful 
there and we do not show the results here. 
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It is useful to understand what part of the flux the dominant contributions come from 
within the log- normal model. The integration kernels are shown in top of figure 1. One can 
see that the mean flux is heavily dominated by the high F region, as expected, caused by 
the presence of r in the denominator of Wp kernel, which leads to divergence for r = 0, 
F = 1. This divergence is cured by the gaussian PDF p(5g\&), but this happens only for 
values extremely close to F = 1, so the preference towards high values of F remains. This 
divergence in itself is not very meaningful: we could also have defined the relevant quantity 
as (1 — F), which would have cured the divergence at F = 1. 

The flux rms a F receives contributions away from the mean flux since W a n kernel 
contains (F— (F)) 2 term. It is again divergent at F = 1 because of r = in the denominator, 
but is also divergent at F = because of F in the denominator of equation 3.13 (this is again 
cured by the gaussian PDF p(5c\cr), but this only happens very close to F = 0). The peak at 
low F is significantly higher than the one at high F. The rms fluctuations are thus heavily 
dominated by the high absorption regions very close to F = 0. 

In contrast to (F) and a F , the bias bp and velocity bias b v are more broadly distributed 
over the entire range < F < 1. The velocity bias b v is simply given by integrating the 
gaussian probability distribution p(5g\&) over dF. The integrand is positive definite over the 
entire range of F. Note that the models that give the same PDF also give the same b v , as 
shown in equation 3.6. We find b v = —0.1 ± 0.01 is the best prediction of this model given 
the observed PDF. 

The bias bp has a similar behaviour for low F, high 5 (or So), where (z^2<5+l)/(l+5) ~ i>2- 
However, for high F, corresponding to 5 < — v% , we see there is a zero crossing of Wb F , so 
those regions contribute positively to bias, while the region 5 > contributes negatively 

to bias. Effects on bp are thus more complicated since the integration kernel is not negative 
definite everywhere: instead, there is a cancellation of positive and negative contributions 
to the integral. Moreover, approaching F = 1, 6 = — 1, the contributions become large 
because of 5 + 1 term in the denominator. The overall result is that bp is susceptible to 
the flux PDF in the voids, which is close to the continuum and so rather poorly determined 
from observations, suggesting the predictions for bp from flux PDF may be less reliable than 
predictions for b v . 

It is instructive to investigate if another nonlinear transform can cure this sensitivity 
to void regions. Since we want to suppress the contribution where F ~ 1 the simplest 
example is to define a new field H = F - F 2 /2, whose kernel is W blI = W bp {l - F) and 
Wb Hv = Wb Fv (l — F). This is also shown in figure 1, and shows that it has achieved the 
desired effect of suppressing the contribution from the void regions. However, it also reduces 
the absolute value of the bias and velocity bias, roughly by 0.05. Such transformations may 
thus be useful if a more robust prediction of bias is needed, although they are likely to 
increase the noise. 

3.2 Effects of scatter 

The above derivation assumes a deterministic relation between the density field and the 
nonlinear transformation. Often the relation is stochastic. For example, galaxies can be 
viewed as a stochastic (Poisson) sampled tracer of a nonlinear transform of the underlying 
density field and the shot noise term is added to this relation, e.g. 5 g = b\5 + &2<5 2 + e. 
Here the scatter is additive and added after the nonlinear transformation. We define e as a 
random variable uncorrelated with the density field, with (e) = and (e 2 ) = a 2 . In Lyman-a 
forest, neutral hydrogen density is a stochastic tracer of the nonlinear transformation of the 
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density field, which we can model as r = r(l + e), with f = (1 + S) a , after which another 
nonlinear transformation, F = exp(— r), takes place. We assumed a multiplicative form of 
scatter to avoid the unphysical situation r < 0, valid as long as a e <C 1. Thus depending 
on the model we can add the scatter either before or after the nonlinear transformation. If 
scatter is uncorrelated with the density field, de/d5i = 0, then there is no contribution from 
the scatter to the calculation of the large scale bias above. 

However, scatter still affects the PDF or the moments of the transformed variable and 
one must correct for this if the moments or the PDF of the transformed variable is used 
as a constraint. For example, in the case of Lyman-a forest we use the observed PDF to 
determine the bias. If we denote F = exp(— f) as the flux in the absence of scatter and 
F = exp(— r) the observed flux with scatter, then since F = F 1+e and dF/dF = F € (l + e) 
we can deconvolve the observed PDF dn/dF to find, 

^/(^L a+ / <(i+£)!,(£K)<i£ - pi5) 

where p(e\a e ) = (2ira 2 )~ 1 ^ 2 exp(— e 2 /2a 2 ) is the PDF for scatter, assumed to be gaussian. 
Simulations suggest a t ~ 0.1 — 0.3 at z = 3 and 0.3-0.5 at z = 2 [4]. The result of this 
analysis for cr t = 0.4 is shown in figure 1. We see that the effects on the PDF are modest. 
Once we have the deconvolved PDF we can proceed as in previous analysis. At z = 2.4 we 
find the effect is to increase \bp\ by about 5% and \b v \ by 3%. For lower values of a € we find 
even smaller effects. 

3.3 Predictions for bp and b v from PDF of Lyman-Q forest 

We now turn to the observational constraints on bp, b v and (3 given the observed PDF. 
The bias and (3 are a function of the two parameters of the model, a and A. We do so by 
evaluating the expressions in equations 3.5-3.6 using observed PDF in [22, 23]. At z = 2.4 
we found that log-normal PDF is a reasonably good fit to the observed one, giving b v = —0.1 
Hence we expect the results to be similar to the best fit PDF in figure 1. This is indeed 
the case: we find b v = -0.09 for [22] PDF and b v = -0.095 for [22] PDF. We note that b v 
is independent of the nonlinear transformation parameters. Bias values depend on the two 
parameters of the transformation. We first compare the prediction to the log-normal model. 
For bias we find, for A = 0.3 and a = 1.6, bp = —0.13 and bp = —0.14 for the two PDFs, 
respectively, compared to bp = —0.11 for log-normal model. It is expected that there will be 
more of a difference between these and the log-normal model since in the void region where 
F ~ 1 the PDF is poorly measured, yet it makes a large contribution with the opposite sign. 
The amplitude A is inversely proportional to UV background photoionization rate T, 

A = 0.96[(1 + z)/4] 4 - 5 T 4 - a7 r:{ 2 , (3.16) 

where we assumed h = 0.7, Cl m = 0.27 and Q^h 2 = 0.0225 and expressed photoionization rate 
r in units of 10~ 12 s~ 4 and temperature T in units of 10 4 iT. Typical values are T_i2 ~ 0.5 — 2 
and T4 ~ 2, relatively independent of redshift between 2 < z < 3 [24], making A to be 
increasing with redshift due to (1 + z) 4,5 dependence. Using equation 3.16 and the expected 
range of UV background amplitude one finds the allowed range is 0.15 < A < 0.6 at z = 2.4. 
The value that fits best the mean flux in hydrodynamic simulations is A = 0.17 [25]. Hence 
A = 0.3 is possibly too high for the observed mean flux. Varying the parameter A while 
fixing a = 1.6 at z = 2.4 we find bp ~ -0.17 at A = 0.15 and bp ~ -0.07 at A = 0.6 
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(since the differences between the two published PDFs is small we simply quote the average 
between the two). Varying a also affects the predicted vaues of bias: for A = 0.3 we find it 
varies from 6p ~ —0.07 at a = 1.2 to bp ~ —0.18 at a = 2.0. 

For A = 0.17 and a = 1.6 we find 6p ~ —0.17 and in combination with b v ~ —0.095 
this gives (3 ~ 0.6. RSD j3 can be as low as 0.5, since bp = 0.2 is at the upper end of 
predicted values, and as high as 1.5, since bp can be as low as -0.06, albeit possibly at an 
unrealistically high A or unrealistically low a. Observed values suggest f3 = 0.8 ± 0.2 [6]. 
Observations constrain best the parameter combination bp + b v and the observed value at 
z = 2.4 is bp + b v = —0.30 ± 0.01 [6] (note that in recent literature it has become standard 
to divide the flux by the mean flux F, while we use the flux itself, hence we multiply the 
published value by F and interpolate to z = 2.4 from z = 2.25 using the measured redshift 
evolution). At A = 0.17, a = 1.6 our prediction is bp + b v = —0.26, slightly below the 
observed value. However, a 10% increase of a = 1.6 or a 10% decrease of A = 0.17 can 
accommodate the observed value. For comparison, the linearized prediction from the term 
linear in 5 is = Aaexp (—A), which, using a = 1.6, gives -0.35 for A = 0.3 (compared to 
our predicted value of -0.13) and -0.5 for A = 0.5 (compared to our prediction of -0.09), while 
for A = 0.17 the linearized model predicts -0.14. We see that the linearized model grossly 
differs from our model for high values of A where the linear approximation is inadequate, 
while for low A the two are in a better agreement. 

The most robust prediction we make is for b v ~ —0.1, since it does not depend on 
any of the transformation parameters, just on the PDF. Comparing to the observations, for 
b F + b v = -0.3 ±0.01 and /3 = 0.8 ±0.2 [6] we find observations suggest b v = -0.13 ±0.03 at 
z = 2.4, which is within one sigma of our prediction. This prediction is robust in the sense 
that it is an analytic prediction of equation 3.6 with no dependence on the transformation 
parameters, only on the flux PDF. At this redshift the published PDFs agree with each other. 
This predicted value is lower than the value measured in simulations of [26], where b v ~ 0.17 
and /3 ~ 1.6. This discrepancy could be due to the simulations not matching the observed 
PDF, or due to additional nonlinear effects that need to be included in our model. 

While the agreement between the predictions and observations is remarkably good, it 
should be pointed out that the real data contain absorbers with Lorentzian wings that cannot 
be modeled as a simple nonlinear transformation F = exp(— ^4(1 + S) a ). These high column 
density systems such as Damped Lyman a systems and Lyman limit systems increase the 
absorption and it is likely that these high density regions increase the bias (in absolute sense) . 
Removing these regions in the data had almost no effect on bp + b v [6], but only a small 
subset of these regions was identifiable in the noisy SDSS data. It remains an open issue 
how much these regions affect the bias. In general, only simulations can address these issues 
in detail. The value of our model is that it identifies the physical effects affecting the bias 
determination, while its quantitative predictions need to be tested in more detail against 
simulations and observations. 

At higher redshifts the published PDFs agree less well with each other, a consequence 
of the fact that the continuum is less well defined. For canonical value of T_i2 = 0.5 — 2 and 
using the two PDFs we predict bp ~ —0.2 ±0.05 at z=3. The predicted values of f3 are similar 
to z = 2.4, around 0.6 to 0.8 in most cases although values outside this range are possible. 
There are no published values for /3 from the data, but extrapolating the measurements of [6] 
to z = 3 one finds bp + b v ~ 0.35, compared to our prediction of 0.35 ± 0.1. This should be 
compared to the linearized prediction of bp ~ —0.35 to -1. The conclusion from this is that 
the analysis gives values of bias and velocity bias in a reasonable agreement with published 
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values, while the simple linearized model grossly overestimates the bias. 

The most important advantage of our model is the possibility to analytically study 
the influence of the physical parameters on the large scale bias. We find that the results 
are rather sensitive to small changes in the two model parameters. This is not surprising 
given the highly nonlinear nature of the transformations. It means that one can determine 
a combination of the two with the large scale bias measurements of bias. Since the two 
parameters also determine the flux PDF a combination of the two measurements provides 
a way to determine both parameters. In particular, assuming A is determined by the mean 
flux then our analysis allows one to measure a from the bias, a complementary method to 
the one used in [27] where the PDF was used to determine a by comparing the simulations 
to the observations. The large scale velocity bias is fully determined by the observed PDF, 
so this allows a robust test of our model predictions. Before we applying this to the data we 
should test our predictions against simulations, so a more detailed analysis of this type will 
be left for the future. 

3.4 Predictions for primordial nongaussianity 

As discussed above primordial nongaussianity induces an additional bias, which can be scale 
dependent for local (equation 2.12), orthogonal or other configurations. We parametrize the 
overall effect in equation 2.19 with the bias 6/ NL . Equation 2.19 applied to our Ly-a forest 
model gives 

dF 2(l-ag) ,, 17 v 

b ^ = 2{6 i6 } = ^^r bF - (3 - 17) 

We see that the prediction for £>/ nl /&f is determined by the values of a and j3. For observed 
value of /3 = 0.8 and for a = 1.6 we find 

^ ~ -0.9. (3.18) 
of 

Note that in the redshift space the observed power constrains better bp + b v [6], so a more 
relevant number to compare is 

0.5. (3.19) 

b F + b v 

Equation 3.18 should be compared to the case of biased galaxies, where that ratio is 
given by = 25 c (b — l)/b ~ 3.37(6 — l)/b, where 5 C = 1.68 (redshift space distortions 
are small for the highly biased objects). The first thing to note is that the sign of the 
effect in Ly-a forest is negative for the chosen parameters, while the effect is positive for 
biased tracers with b > 1: primordial nongaussianity with positive /nl supresses power in 
Ly-a forest. We also see that the relative effect of primordial nongaussianity in Ly-a forest 
is significantly reduced relative to the highly biased tracers, since 0.5 or 0.9 is much less 
than 3.37. This prediction in fact sensitively depends on the parameters a and /3, as seen 
in equation 3.17. The primordial nongaussianity effect vanishes if j3 = a -1 , since this is 
the same as the linearized limit F hn = Aa exp (— A)5, where the primordial nongaussianity 
signatures do not appear in 2-point correlations. The effect changes sign for (3 < a" 1 relative 
to the sign above. On the other hand, it should be pointed out that the overall effect scales 
by the inverse of the growth rate D(z) (equation 2.13), which for Ly-a forest redshifts can 
be a factor of 2-3 larger relative to z = 0. 

Since the predictions for Ly-a forest are sensitive to the parameters and very close to 
zero it is worth exploring if a further nonlinear transformation of Ly-a forest can change 
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the nongaussian bias relative to gaussian bias. For example, for H = F — F 2 /2 one finds 
bn = —0.07, b v = —0.04 and 6/ NL = —0.04, so that bf NL u/bH = 0.6, reversing the sign of the 
effect. Similarly, defining a new field as G = (F — F) 2 gives bo = 2[(1— F)bp — bu], so that for 
F = 0.8 one has ^/ nl g/^g = 1-2, which is also positive in sign and a factor of 2 larger. If we 
want a stronger contrast we can arrange the bias to vanish. For example, for K = 0.6-F — H 
we find bf Nh K/^K = oo, because the large scale bias bx vanishes for this combination. We can 
also choose a different combination such that (3 or 6(1 + j3) vanishes (but we cannot make all 
the large scale power in redshift space vanish since we cannot simultaneously vanish density 
and velocity bias). Still, such transformations may boost the primordial nongaussianity signal 
relative to the gaussian signal on large scales. This is related to the multi-tracer method of 
canceling sampling variance [20], where two tracers with different bias values are combined 
in a way that the large scale fluctuations in the density field cancel: this in fact happens if 
one considers the linear combinations of the two tracers where the large scale bias vanishes, 
in which case one is left with the large scale primordial nongaussianity bias, which does not 
vanish. Such nonlinear transforms do not necessarily enhance the signal itself and typically 
increase the noise properties of the new field, so only in the sampling variance limit one gains 
using the sample variance canceling techniques. We do not expect to be in such a limit for 
Ly-a forest with current data sets [28, 29]. 

Yet another potential advantage of these nonlinear transformations is to test against 
systematics. On large scales there are numerous sources of additional power in Ly-a forest 
such as UV background fluctuations or He reionization signatures [30, 31], which could mimic 
the scale dependent bias effects of primordial nongaussianity. However, primordial nongaus- 
sianity changes in a predictable way under the nonlinear transformation. For example, if 
one observes an unexplained large scale power in Ly-a forest compatible with primordial 
nongaussianity one can devise a nonlinear transformation where primordial nongaussianity 
vanishes, while one would not expect the rest of the effects to vanish. A more detailed analysis 
to see how much can be gained by such methods is beyond the scope of this paper. 

4 Conclusions 

In this paper we discuss the large scale clustering properties of nonlinear transformations of 
density field. We show that on large scales these can be viewed as biased versions of the 
density field itself. We present an analytic derivation of bias, velocity bias and primordial 
nongaussianity bias. The interesting aspect of our calculation is that the large scale bias can 
be expressed entirely in terms of the final PDF and the physical parameters of the nonlinear 
transformation, allowing one to determine it from the observed PDF without the need for 
simulations. The resulting bias depends on the nature of the transformation and with a 
suitable nonlinear transformation one can design a field with very different large scale bias 
than the original field, including bias of zero. 

The primary application in this paper is to the ionizing equilibrium model of Lyman-a 
forest: we derive the bias, velocity bias and primordial nongaussianity bias of Ly-a forest 
flux as a function of redshift from the observed PDF. Velocity bias has a very simple ex- 
pression (equation 3.6) that depends only on the observed PDF and its predicted value is in 
a good agreement with observations of [6] . We derive the bias as a function of parameters 
of transformation, which are dominated by the UV background amplitude and the slope of 
the temperature-density relation. Assuming their fiducial values the values of bias and fi are 
in a reasonable agreement with observations of [6]. We find that the bias has a significant 
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contribution from the void regions, where lognormal PDF differs from the measured one, as 
seen in figure 1 . Another nonlinear transformation, e.g. to H = F — F 2 /2, can reduce 
the sensitivity to voids. We find that the primordial nongaussianity bias for the observed 
f3 ~ 0.8 has the opposite sign than highly biased galaxy tracers, but is relatively small. A 
small change in the parameters can change its sign or make it zero. 

There are several generalizations and applications of the model that can be pursued. 
First, we used the linear growth evolution to derive the coupling between long and short 
wavelength modes. This may break down when a 2 = (5 2 ) 3> 1. It would be useful to extract 
from simulations the general value of coefficients d5 n /d5i as a function of scale, and use these 
in the calculations above. With this one should be able to fully compute the bias for any 
nonlinear transformation and connect it to its PDF. One can also extend the calculation 
to the next order and investigate second order bias, 62 = (d 2 T / d5f) /2. For example, one 
could design a transformation where second order bias vanishes, therefore suppressing the 
scale dependence of the bias. Similarly, one can also design a transformation where velocity 
bias vanishes, making the redshift space correlation function isotropic on large scales. Or, 
one could design a transformation where the bias vanishes, in which case all the large scale 
correlations would come from velocities. 

Having established that a nonlinear transformation traces the large scale modes, but 
with a different bias, another interesting question worth further investigation is whether such 
transformations can increase the signal to noise (S/N) of the power spectrum measurement. 
On large scales this is proportional to b 2 /a 2 , where a 2 is the variance of the field, includ- 
ing observational noise. While an increase in bias increases the signal, the accompanying 
nonlinear transformation may also increase the noise, so it is not clear if S/N can increase. 
For example, taking a log of flux F gives optical depth r = — In F, which for a = 2 has 
bs T ~ 2 — 3 (equation 2.10), compared to bp ~ —0.15, a huge increase in bias. However, 
the corresponding nonlinear transformation also greatly enhances the noise, both detector 
and photon noise, specially for F ~ 0, as well as Poisson like noise caused by high neutral 
hydrogen column density sources in the optical depth, which makes this example likely to be 
impractical. There may however be other nonlinear transformations that perform better. A 
detailed analysis is beyond the scope of this paper. 

Another application of the method is to multi-tracer methods: in [19, 20] it was argued 
that combining two tracers with different bias parameters can reduce the sampling variance 
on certain quantities, such as the RSD (3 or nongaussianity parameter f n [. This is because 
if the two tracers both trace the underlying long wavelength 5i, taking the ratio of the two 
eliminates 5i itself, and so eliminates the dominant source of error on large scales, which is 
the stochastic nature of 5i (i.e., sampling variance). At the same time this ratio contains 
some useful cosmological information, such as sensitivity to primordial nongaussianity [20]. 
Here we have shown that the second tracer with a different bias can be simply obtained by 
a nonlinear transformation of the field itself and that takig a linear combination of the two 
can result in vanishing of the large scale bias, effectively achieving the sampling variance 
cancellation. These methods are of course simply reduced versions of combining N-point 
statistics, but are particularly simple to understand and analyze. In this view the nonlinear 
transformation converts the higher order correlations of the density field into the 2-point 
correlations of the nonlinear transform. For example, if we square the density field and 
correlate with the density field itself we obtain a reduced version of a 3-point function, which 
we can then compare to the 2-point function of the density field. Both will trace the long 
wavelength modes on large scales, but the first one will be biased relative to the density 
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field itself. Comparing the two thus eliminates the sampling variance, while preserving some 
cosmological information in the bias itself: for example, its sensitivity to primordial local 
nongaussianity has been shown in e.g. [15]. While this example involves comparing 2 and 
3 point functions, other combinations may be more effective in extracting the information 
optimally. 

Testing the results against systematics may be another application of such nonlinear 
transforms. This is specially true for primordial nongaussianity, which enhances the power 
on very large scales. Other effects may achieve the same effect: for Ly-a forest the additional 
enhancement of power could come from UV background fluctuations or He reionization, for 
galaxies the large scale power could be enhanced due to the effects from our own galaxy, 
such as star-galaxy separation or extinction. In most cases it may be difficult to separate 
these effects from the primordial nongaussianity. However, one can use a nonlinear tranform 
to change the primordial nongaussianity bias and using the expressions derived here we can 
predict what their large scale bias is. For example, as discussed in this paper one can devise 
a noninear transform where primordial nongaussianity bias vanishes, while one would not 
expect the rest of the effects to vanish, so one can separate the two effects. 

In summary, the results of this paper can be used to develop an analytic understanding 
of the large scale bias of nonlinear transformations and its sensitivity to parameters of the 
transformation. This is required if we want to achieve a better theoretical understanding of 
bias of Ly-a forest and other nonlinear tracers of matter density. Our results allow one to 
design nonlinear transformations with a nearly arbitrary large scale bias, which can be used 
to optimize the extraction of cosmological information from the cosmological observations. 
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